function VARbc = doCholVARbootstrap_CS(VAR,nboot,clevel )

res      = detrend(VAR.res,'constant');
for jj=1:nboot;

         Tmax = max(VAR.Ti);                                             
         cal  = (min(VAR.Year):1:max(VAR.Year) )' ;
        rr    = 1-2*(rand(Tmax,1)>0.5);                                    
        resb = zeros(VAR.T,VAR.Ny);
        for nn =1:VAR.T;
            resb(nn,:) = rr(find(cal==VAR.Year(nn))) .* res(nn,:) ; 
        end
    
           t  = 1;
           Yb = zeros(VAR.T,VAR.Ny);
           Xb = zeros(VAR.T,VAR.Ny*VAR.p); 
           
        for i=1:max(VAR.ii);                                                
            Xb(t,:) = VAR.X(t,1:VAR.Ny*VAR.p) ;                                
            for j=t:(t+VAR.Ti(i)-1);                                 
                 Yb(j,:)   = Xb(j,:)*VAR.bet(1:VAR.Ny*VAR.p,:) + resb(j,:);           
                 if j<t+VAR.Ti(i)-1;
                 Xb(j+1,:) = [Yb(j,:), Xb(j,1:((VAR.p-1)*VAR.Ny))];  
                 end
            end 
            t = t + VAR.Ti(i);
        end
        
        % run VAR with this boostrappred data
        VARBS.p=VAR.p; VARBS.Ny=VAR.Ny; VARBS.N=VAR.N; VARBS.irhor=VAR.irhor; VARBS.Perm=VAR.Perm;
        VARBS.normalize=VAR.normalize;
        VARBS.Y   = Yb;
        VARBS.X   = Xb ;                             
        VARBS.bet = [VARBS.X ones(VAR.T,1)] \ VARBS.Y;              
        VARBS.res = VAR.Y - [ VAR.X ones(VAR.T,1)]*VARBS.bet;
      
        VARBS.Sigma = VARBS.res'*VARBS.res / ( VAR.T - VAR.Ny*VAR.p - 1);   
        % compute Cholesky impulse response function
        VARBS   = irfvar(VARBS);
        
        VARbs.irf(:,:,:,jj)   = VARBS.cholIRF;
        clear VARBS
end 

% IRF confidence intervals
VARbc.CholirsH1=quantile(VARbs.irf,(1-clevel(1)/100)/2,4);
VARbc.CholirsL1=quantile(VARbs.irf,1-(1-clevel(1)/100)/2,4);
VARbc.CholirsH2=quantile(VARbs.irf,(1-clevel(2)/100)/2,4);
VARbc.CholirsL2=quantile(VARbs.irf,1-(1-clevel(2)/100)/2,4);
VARbc.Cholirsmean=mean(VARbs.irf,4);
VARbc.Xb     =Xb;
% IRF cumulative confidence intervals
VARbc.CholirsCumH1=quantile(cumsum(VARbs.irf,1),(1-clevel(1)/100)/2,4);
VARbc.CholirsCumL1=quantile(cumsum(VARbs.irf),1-(1-clevel(1)/100)/2,4);
VARbc.CholirsCumH2=quantile(cumsum(VARbs.irf),(1-clevel(2)/100)/2,4);
VARbc.CholirsCumL2=quantile(cumsum(VARbs.irf),1-(1-clevel(2)/100)/2,4);

end

